# Analyzing metabolomics data for Mutographs ----
# 01_data_management:
# Producing the tables that will be used in the rest of the analysis

rm(list = ls())

# 0 - Definition of libraries, paths and functions ----

Sys.setlocale("LC_TIME", "C")

library(stringr)
library(tidyr)
library(Hmisc)
library(MASS)
library(data.table)
library(dplyr)
library(forcats)
library(tibble)
library(openxlsx)
library(ggplot2)
library(cowplot)
library(scales)
library(patchwork)
library(ggpubr)

# 1 - Data import and management ---- 
metabo_RAW <- read.csv("data/metabolomics_normalized_data.csv", check.names=FALSE)
features <- names(metabo_RAW)[4:2395]

# Sparseness of data
# Checking the number of missing (NA) values per feature
sparseness <- colSums(is.na(metabo_RAW)) %>% 
  as.data.frame() %>% rename(N0 = ".") %>% 
  rownames_to_column() %>% 
  filter(!rowname %in% c("donor_id", "batch", "acq_order")) %>% 
  arrange(desc(N0)) %>% 
  mutate(P0 = 100*N0/901, 
         category = cut(P0, breaks = c(min(P0)-1,5,25,50,75,95,max(P0)+1), labels=c("<5%","5-25%","25-50%","50-75%","75-95%",">95%")))

# Sorting features by completeness
features_sorted <- arrange(sparseness, N0)$rowname 

# Pre-processing the data
# Removing the features that are appearing only in one batch and not the other one (possibility of a technical glitch)
rm_features_batch <- metabo_RAW %>% 
  group_by(batch) %>% 
  summarize(across(all_of(features), ~sum(!is.na(.)))) %>% 
  t() %>% 
  as.data.frame() %>% 
  rownames_to_column() %>% 
  filter(V1 == 0 | V2 == 0) %>% 
  pull(rowname)
# Only 21 features wi

sel_features <- features_sorted[!features_sorted %in% rm_features_batch]

# Calculating correlation between features
metabo_pearson <- cor(select(metabo_RAW, all_of(sel_features)), use = "pairwise.complete.obs")
saveRDS(metabo_pearson, "output/metabo_pearson")

# Recursive filtering elimination of features 
# In a cluster of features (r > 0.85), only selecting the less sparse feature
# NB: Can remove the printing lines to add efficiency to the code
list_excluded <- c()
matrix_corr <- metabo_pearson
proxy_table <- data.frame()
for (metabolite in sel_features){
  if (!metabolite %in% list_excluded){
    print(str_c("Evaluating ", metabolite))
    df <- matrix_corr %>% 
      as.data.frame() %>% 
      rownames_to_column() %>% 
      mutate(N = row_number()) %>% 
      filter(abs(!!sym(metabolite)) > 0.85 & rowname != metabolite) %>% 
      select(rowname, all_of(metabolite), N) %>% 
      mutate(feature = metabolite) %>% 
      rename(removed_feature = rowname, pearson = metabolite)
    
    excluded <- df$removed_feature
    idx_excluded <- df$N
    if (NROW(df != 0)){
      matrix_corr <- matrix_corr[-idx_excluded,-idx_excluded]
      list_excluded <- c(list_excluded, excluded)
      print(str_c("Excluding ", paste(excluded, collapse=", ")))
      print(str_c("Number of remaining features: ", NROW(matrix_corr)))
      proxy_table <- bind_rows(proxy_table, select(df, feature, removed_feature, pearson))
    } else {print("No features excluded")}
  } else {print(str_c(metabolite, " was already excluded"))}
}
rm(df, excluded, idx_excluded)
# In proxy_table: listing all metabolites - the one kept in the analysis and the one removed - with the corresponding Pearson's r
saveRDS(proxy_table, "output/proxy_table")
saveRDS(list_excluded, "output/list_excluded")

# List of final features to be included in the analysis
metabolites <- sel_features[!sel_features %in% list_excluded]

# Features are log-transformed and Pareto scaled (mean-centered and divided by the square root of the standard deviation of each variable)
metabo_PROC <- metabo_RAW %>% 
  # Removing one-batch features and highly correlated features
  select(-all_of(c(rm_features_batch, list_excluded))) %>% 
  # All metabolites to numeric
  mutate(across(all_of(metabolites), as.numeric)) %>% 
  # Negative to NA
  mutate(across(all_of(metabolites), ~ifelse(.<0,NA,.))) %>%
  # Missing to 1/5 of minimum
  mutate(across(all_of(metabolites), ~replace_na(.,1/5*min(.,na.rm=T)))) %>%
  # Log transformation
  mutate(across(all_of(metabolites), ~log10(.))) %>% 
  # Pareto scaling
  mutate(across(all_of(metabolites), ~(. - mean(.,na.rm=T))/sqrt(sd(.,na.rm=T))))

# Adding metadata and signatures
# Epidemiological and clinical variables
rcc_CORE <- read.csv("data/meta/MutWP1_RCC_core_data_Manuscript_v3_1.csv", stringsAsFactors = F) %>% select(-lastupdt, -log_updt)
rcc_TOB <- read.csv("data/meta/MutWP1_RCC_hypothesis_specific_tobacco_Manuscript_v3_1.csv", stringsAsFactors = F) %>% select(-lastupdt, -log_updt)
rcc_ALC <- read.csv("data/meta/MutWP1_RCC_hypothesis_specific_alcohol_Manuscript_v3_1.csv", stringsAsFactors = F) %>% select(-lastupdt, -log_updt)
rcc_TUM <- read.csv("data/meta/MutWP1_RCC_tumour_specific_Manuscript_v3_1.csv", stringsAsFactors = F) %>% select(-lastupdt, -log_updt)

# Signature attribution - De novo SBS, COSMIC SBS, COSMIC DBS, COsMIC ID, total mutational burden
sig_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_denovo_SBS1536_abs_mutations.csv") %>% rename(donor_id=X)
cosmic_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_COSMIC_SBS96_abs_mutations.csv") %>% 
  rename(donor_id=X, SBS1536A_cosmic = SBS1536A, SBS1536B_cosmic = SBS1536B, SBS1536F_cosmic = SBS1536F, SBS1536I_cosmic = SBS1536I)
dbs_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_COSMIC_DBS78_abs_mutations.csv") %>% rename(donor_id=X)
id_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_COSMIC_ID83_abs_mutations.csv") %>% rename(donor_id=X)
total_mutburden <- read.csv("data/sigs/output_RCC_Manuscript_denovo_SBS1536_stat_info.csv") %>% 
  rename(donor_id=Sample, mutburden=Mutational.burden) %>% select(donor_id, mutburden)

all_RCC <- rcc_CORE %>% left_join(rcc_TOB) %>% left_join(rcc_ALC) %>% left_join(rcc_TUM) %>%  
  left_join(sig_attrib) %>% left_join(total_mutburden) %>%
  left_join(cosmic_attrib) %>% left_join(dbs_attrib) %>% left_join(id_attrib)

sigs <- names(sig_attrib)[2:13]
sigs_cosmic <- names(cosmic_attrib)[2:15]
dbs_cosmic <- names(dbs_attrib)[2:6]
id_cosmic <- names(id_attrib)[2:10]

# Creating dataset with features + metadata + signatures
metabo_MAN <- metabo_PROC %>% 
  left_join(all_RCC) %>% 
  filter(donor_id != "PD47592a") %>%  # Removing outlier from further tests
  mutate(across(where(is.character), ~na_if(.,"Missing"))) %>% 
  mutate(country = relevel(factor(country), ref = "Czech Republic"),
         sex = relevel(factor(sex), ref="Male"),
         bmi = as.numeric(bmi),
         batch = relevel(factor(batch), ref = "1")) %>%
  mutate(across(
    any_of(c(sigs,sigs_cosmic,dbs_cosmic,id_cosmic,"mutburden")),
    .fns = list(cat = ~cut2(.,g=2), int = ~qnorm((rank(.,na.last="keep")-0.5)/sum(!is.na(.))), logdelta = ~log2(.+1)),
    .names = "{.col}_{.fn}"))

saveRDS(metabo_PROC, "output/metabo_PROC")
saveRDS(metabo_MAN, "output/metabo_MAN")